// WindowAnalysis.cpp: implementation of the WindowAnalysis class.
//
//////////////////////////////////////////////////////////////////////

//#include "stdafx.h"
#include "WindowAnalysis.h"
#include<iostream>
#include<fstream>
#include<complex>
using namespace std;

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

WindowAnalysis::WindowAnalysis()
{
	expval = NULL;
	 s	   = NULL;
     x	   = NULL;
	 framecount=0;


}

WindowAnalysis::~WindowAnalysis()
{
	
	delete []expval;	
	delete []s	  ;
	delete []x	  ;
	//printf("Destructor of WindowAnalysis ");
	//fflush(stdout);
  //delete []CRegion;


}

WindowAnalysis::WindowAnalysis(int len , double thold ,long intron , long exon)
{

	lenWindow=len;
	thersHold=thold;
	exonCutOff=exon;
	intronCutOff =intron;
	
	s	   = NULL;
    x	   = NULL;

	framecount=0;
//	cout<<"ThresHold Value "<<thersHold;
//	cout<<"Intron :"<<intronCutOff<<" Exon :"<<exonCutOff<<endl;

    expval=new complex<double>[lenWindow];

}


void WindowAnalysis::settherHold(double thold)
{
	thersHold =thold;
//	 cout<< "ThresHold Value "<<thersHold;

}

void WindowAnalysis::setWindowLen(int len)
{
    lenWindow=len;
	delete [] expval;
	expval=new complex<double>[lenWindow];

}
void  WindowAnalysis::setIntronExon(long intron , long exon )
{
	 exonCutOff=exon;
	 intronCutOff =intron;


}
int WindowAnalysis::doFt(const char  * strSeq , char * logfile)
{


const double F3 = 1.000/3;
const double PI = 3.14159265;


//ofstream	out("output" );
ofstream	out(logfile ,ios::app);

long lenStr = strlen(strSeq);

if(lenStr<lenWindow)
{
	out<<"Sequence length is less than windowlength."<<endl;
	return 1;
}

cout<<"\n Sequence Length :"<<lenStr<<endl;


int result= findexp(lenWindow);

//long *CRegion = new long[lenStr-lenWindow+ 1];
//long *s	    = new long[lenStr-lenWindow+ 1];
//double *x	= new double[lenStr-lenWindow+1];

delete [] s;
delete [] x;
s	    = new long[lenStr-lenWindow+ 1];
x	= new double[lenStr-lenWindow+1];


char nuclt[5]= "ATGC";

double winsq		= (1.0/lenWindow)+(1.0/(lenWindow*lenWindow));
long  wincube	    = lenWindow*lenWindow*lenWindow;


int countA =0;
int countC =0;
int countG =0;
int countT= 0;

complex<double> s1(0.0,0.0);

double sq=0.0;
double s2=0.0;
double sbar =0.0;
double p=0.0;

framecount=0;



for(long k = 0 ; k < lenStr-lenWindow+1 ; k++ )

{
	sq		=0.0;
	s2		=0.0;
	sbar	=0.0;
	p		=0.0;

	// Frame wise  count of A C G T seperately


if(k==0)
{
	for(long i =k; i < k+lenWindow; i++ )
	{

	if( strSeq[i] == 'A') 
		{ countA++; continue;}
  	if( strSeq[i] == 'G') 
		{ countG++; continue;}
  	if( strSeq[i] == 'C') 
		{ countC++; continue;}
  	if( strSeq[i] == 'T') 
		{ countT++; continue;}

	} //end of count for look.

}else
{

switch(strSeq[k-1])
{
case 'A' :
    { countA--; break;	}
case 'T' :
    { countT--; break;	}
case 'G' :
    { countG--; break;	}
case 'C' :
    { countC--; break;	}
}

switch(strSeq[k+lenWindow-1])
{
case 'A' :
    { countA++; break;	}
case 'T' :
    { countT++; break;	}
case 'G' :
    { countG++; break;	}
case 'C' :
    { countC++; break;	}
}

 
}

//out<<"Count A  :"<<countA <<" Frame count :"<<k<<endl;
//out<<"Count T  :"<<countT <<" Frame count :"<<k<<endl;
//out<<"Count G  :"<<countG <<" Frame count :"<<k<<endl;
//out<<"Count C  :"<<countC <<" Frame count :"<<k<<endl;

 sq = 0.0+ pow(countA,2)+pow(countG,2) + pow(countT,2) + pow(countC,2);

// for each Nucleotide. A T G C


	for(int i = 0; i <4 ; i ++ )
	{
		s1= complex<double>(0.0,0.0);


		for (int bi=0; bi<lenWindow; bi++)
		{
				if(strSeq[k+bi] == nuclt[i] )
					s1= s1+expval[bi];
		}

		s2 = s2+ (abs(s1)*abs(s1));

	}


//	to calculater power peak 

		sq   = sq/wincube;
		sbar =(winsq-sq);
		sq   = sbar * (lenWindow/2.0);

// to calcuater signal to noise ratio

		s2 = s2/(sq*lenWindow*lenWindow);

		p= s2/sbar;


	s[framecount]=k+(lenWindow/2);
	x[framecount]=p;

	out<<s[framecount]<<"                "<<x[framecount]<<"\n";
    framecount++;


} // End of Outer for loop(Frame wise)

//return(0);
//} 
// end of Ft function..
//void WindowAnalysis::findCodingRegion(long seqfileId , const char *  seqId)
//{

long *CRegion = new long[framecount];
//memset(CRegion , 0 ,(framecount+1)*sizeof(long));

int    flag=0; 
long i =0 , m =0;

	for (m =0; m<framecount; m++)
	{
       
		if(x[m]>=thersHold)
		{
			if(flag == 0)
			{
				
				CRegion[i] =s[m];
				//out<< "\n START :" <<s[m]<<" : "<<x[m];
				//i++;
				flag =1;
			}

		}

        else
        {
            if(flag==1)
       	     {
               if(x[m]<thersHold)
               {
				 	
					 	 out<< "\n START :" <<CRegion[i];
						 out<<"\t\tEND: "<<s[m-1]<<" : "<<x[m-1];		 
						i++;
						CRegion[i]=s[m-1];
						i++;
				}
				flag=0;

			 }
			
        } // end of if.
	}// End Coding region loop. 



if (flag==1)
{
	out<< "\n START :" <<CRegion[i];
	out<<"\t\tEND: "<<s[m-1]<<" : "<<x[m-1];		 
	i++;
	CRegion[i]=s[m-1];
	i++;
}

//delete []CRegion;
// that data is deleted 07/02/02 To inprove performance.

//delete []s	  ;
//delete []x	  ;
//s=NULL;
//x=NULL;
 
out<<"\n Final Coding Regions   :"<<endl;

flag =0;
int code =0;
long sta=0 , stp =0;
structCR obj;
obj.type=' ';
					
for(int j=0 ; j<i ; j +=2)
   {
	if(flag ==0)
	{
		sta = CRegion[j];
		stp = CRegion[j+1];
        if ( (stp-sta)>= exonCutOff ) code =1; else code =0; 	
		flag = 1 ;
    }
	else
	{                      
		if( (CRegion[j]-CRegion[j-1]) >=intronCutOff  )
		{
			if(code ==1)
			{	
				out <<"\nStart :"<<sta <<"  End  :" << stp;  			
				// inset in database  start new set.
				    obj.start=sta;
					obj.end=stp;
					vecCRegion.push_back(obj);
				//insetcodes( sta , stp,seqfileId,seqId);
				sta = CRegion[j];
				stp = CRegion[j+1];
			
				if ( (stp-sta)>= exonCutOff ) code =1; else code =0; 	
				flag = 1 ;
			}
			else
			{ 
				//cout << forget the old set.
				  if ((stp-sta)>= exonCutOff ) 
				  { out <<"\nStart :"<<sta <<"  End  :" << stp;
					obj.start=sta;
					obj.end=stp;
					vecCRegion.push_back(obj);
					//insetcodes( sta , stp,seqfileId,seqId);
				  }
				// may need to be  commented later.
				sta = CRegion[j];
				stp = CRegion[j+1];
				if ( (stp-sta)>= exonCutOff ) code =1; else code =0; 	
				flag = 1 ;
			}
		}
		else
		{

        	if( ( CRegion[j+1]-CRegion[j] ) >= exonCutOff ) code =1;
			stp = CRegion[j+1];
        	flag =1;

		}
	}
 }

if ((code == 1) || ((stp-sta) >exonCutOff) )
{

	out <<"\nStart :"<<sta <<"  End  :" << stp;
	obj.start=sta;
	obj.end=stp;
	vecCRegion.push_back(obj);
	//mapCRegion[sta] =stp;
 	//insetcodes( sta , stp,seqfileId,seqId);		
}


delete []CRegion; 
CRegion=NULL;


vector<structCR>::iterator it;

PowerSpectrum powerObj;

char * strtemp=NULL;

for(it =vecCRegion.begin(); it !=vecCRegion.end(); it++)
{
	//long start =it->start;
	//long end =it->end;

	strtemp = new char[( (it->end - it->start) +1)];

	memset(strtemp , '\0' ,(it->end - it->start) +1 );

	strncpy( strtemp , &strSeq[it->start] ,(it->end - it->start) );

	//strncpy( strtemp , &strSeq[start] , end -start);

    printf("Length : %d \n ", strlen(strtemp));

	printf("start %ld , End  %d  length %d\n " , it->start ,it->end , (it->end - it->start));
	fflush(stdout);

	powerObj.setSequence(strtemp);
	int peaks =powerObj.doPowerSpectrumAnalysis();
	if(peaks==0)
	{
		it=vecCRegion.erase(it);
		if( it !=vecCRegion.end() )
		--it;
		else
			break;
	}
	else if(peaks==1)
	{
		it->type ='S';
	}
	else if(peaks==2)
	{
		it->type ='M';
	}
delete []strtemp;
strtemp = NULL;
}



  //store the coding regions in an array


//int elsize;
 //elsize=vecCRegion.size();



delete []strtemp;
strtemp = NULL;


//delete []CRegion; 
//CRegion=NULL;

//delete []s	  ;
//delete []x	  ;

return(0); 

} // end of Coding region finding.
 

// to find exponential positional values to window size .

int WindowAnalysis::findexp(int WindowSize)
{


complex<double>  i1(0.0,1.0);
const double PI = 3.14159274;
const double F3 = 1.000/3;

for(int k=0; k<WindowSize; k++)
{ 
	
	expval[k] = exp(-i1*2.0*PI*F3*(k+1.0));		

}

return(0);

}

// ///////////////// End //////////////
